Thermal expansion in single-walled carbon nanotubes and graphene: nonequilibrium 

Green's function approach 
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The nonequilibrium Green's lunction method is applied to investigate the coefficient of thermal 
expansion (CTE) in single- walled carbon nanotubes (SWCNT) and graphene. It is found that atoms 
expand about 1% from equilibrium positions even at T = K, resulting from the interplay between 
quantum zero-point motion and nonlinear interaction. The CTE in SWCNT of different sizes is 
studied and analyzed in terms of the competition between various vibration modes. As a result 
of this competition, the axial CTE is positive in the whole temperature range, while the radial 
CTE is negative at low temperatures. In graphene, the CTE is very sensitive to the substrate. 
Without substrate, CTE has large negative region at low temperatures and very small value at high 
temperature limit, and the value of CTE at 300 K is —6 x 10~^ which is very close to a recent 
experimental result, -7 x lO"*^ K"^ [Nat. Nanotechnol. 10, 1038 (2009)]. A very weak substrate 
interaction (about 0.06% of the in-plane interaction) can largely reduce the negative CTE region 
and greatly enhance the value of CTE. If the substrate interaction is strong enough, the CTE will be 
positive in whole temperature range and the saturate value at high temperatures reaches 2.0 x 10~^ 

PACS numbers: 65.80.-|-n, 05.70.-a, 61.46.-w, 62.23.Kn 
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I. INTRODUCTION 



Single-walled carbon nanotubes (SWCNT) and 
graphene are two kinds of novel carbon based mate- 
rials with lots of intriguing electronic and mechanical 



properties 



They have many potential appli- 



cations in nano-devices, where the thermal property 
plays a big role. The coefficient of thermal expansion 
(CTE) is one of the most important nonlinear thermal 
properties. Very recently, the CTE of graphene is 
measured to be -7 x 10"^ at T = 300 Also 

the CTE of multi-walled carbon nanotubes^'^'i^ and 
nanotube bundle o^^'^^i^^'^^ are measured in experiments. 
Theoretically, there are several approaches to calculate 
CTE such as molecular dynamics r^ii^iiiii^ lattice dy- 
namics calculationsr^ii^'^i2ii^ the molecular structural 
mechanics method^^, analytical method^'*, and ab initio 
density-functional theory and density-functional pertu- 
bation theory calculationsi^i^ However, discrepancies 
exist in different approaches. In Ref. [TUlil, the 
CTE is predicted to be negative in low temperature 
region while positive at high temperatures. On the 
contrary, Li and Chou obtain positive value for CTE 
in whole temperature rangci^ In Ref. [2^, CTE for 
graphite is negative at low temperatures and positive at 
high temperatures, and CTE for graphene is negative 
in whole temperature range. In Refs. [T8ll2^ . the axial 
CTE increases from negative to positive with increasing 
temperature; while the minimum value of axial CTE is 
quite different: -1 x lO"'^ K'^ or -5 x 10"^ K'^ and 
this value is achieved at different temperatures: 300 K 
or 500 K in Ref. [H and Ref. [H respectively. In the 
high temperature region, the value for CTE is also quite 



different from one to another, typically about 2.0 x 10 ^ 
K-i in Refs. MM, 4.0 x 10"^ K"! in Ref. [13 or more 
than 5 X lO"'^ K'^ in Ref. [H. In Ref. [H the CTE shows 
only a slight dependence on the length or diameter of the 
SWCNT. While Ref. [H predicts the obvious decrease 
of CTE with increasing diameter. To the best of our 
knowledge, the limited existing theoretical works are 
diverse. 

In this paper, we investigate the CTE in SWCNT and 
graphene sheets by the nonequilibrium Green's function 
(GF) approach, which includes contributions from all 
phonon modes and takes into account the quantum ef- 
fect. Our study shows that even at T = K the aver- 
aged position of each atom deviates about 1% from its 
equilibrium position. This is the result of the combined 
contribution from quantum zero-point motion and non- 
linear interaction. For the CTE in SWCNT, we study the 
competition between lateral bending vibration (negative 
effect), the radial breathing mode (positive effect), and 
the longitudinal optical phonon mode (positive effect). 
As a result, the axial CTE is positive in the whole tem- 
perature range; while the radial CTE has negative value 
at low temperatures and the radial CTE is in general 
one order of magnitude smaller than the axial CTE at 
high temperatures. Our results show that the axial CTE 
is not very sensitive to the size of SWCNT. Yet for ra- 
dial CTE, the thinner SWCNT (with higher aspect ratio 
length/diameter, L/d) has more tense bending vibration, 
thus has smaller value. For graphene, we demonstrate 
that the value of CTE is very sensitive to the interac- 
tion 7 between the substrate and graphene. If there is 
no substrate, graphene has negative CTE at low tem- 
peratures and generally the CTE is very small at high 



temperatures, and the value of CTE at T = 300 K from 
our calculation is a = —6 x 10^® K^^, which agrees with 
the experimental one. When the interaction is nonzero, 
even a very small value (7 — 0.001 eV/ (A^u)) will largely 
shrink the negative CTE region, and greatly enhance the 
value of CTE at high temperatures. If 7 is larger than 
0.1 eV/(A^u), the CTE is positive in the whole temper- 
ature range and saturates at value of 2.0 x 10~^ at 
high temperature. 

The paper is organized as follows. In Sec. II, we present 
the detailed derivations in the GF method. Sec. Ill is 
devoted to the main results and relevant discussions. The 
conclusion is in Sec. IV. 



II. GREEN'S FUNCTION METHOD 
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The definition and notations of GF used in this paper 
are from Rcf. 27,28,29, where relationships between dif- 
ferent versions of GF can also be found. In the present 
investigation of CTE, two types of GF are required : 

G,(r) = -^(T^ufir)), (1) 

GMr,r') = -J(r.<(r)«f(r')), (2) 

where {t) is the vibrational displacement of atom 
j in Heisenberg picture, multiplied by the square root 
of mass of the atom. Tj- is the contour-order opera- 
tor on the Keldysh contour and r is the time on the 
contour. The contour-order operator in Gj{T) is nec- 
essary in the further perturbation expansion of it. In 
the harmonic system without nonlinear interaction, all 
atoms vibrate around the equilibrium position linearly, 
so G^(t) = and the two-point harmonic contour or- 
dered GF G'°j,(t, r') can be exactly calculated: 



G?.(r,r') = -l(T.^f(r)^f(r'))o, 



(3) 



where the subscript indicates a harmonic system. 

These two GFs are defined in the Heisenberg picture, 
however the use of interaction picture is more advanta- 
geous in the study of systems with nonlinear interaction. 
The Gj (r) in the interaction picture can be obtained by 
using the general picture transformations: 



(4) 



where is the vibrational displacement in the interac- 
tion picture. S is the scattering matrix. The nonlinear 
interaction used in the present paper takes the form: 



Hn = —^UlUmUn 

lm7i 



E'^opqr 
4 

opqr 



UoUpUqUr, (5) 



where kimn and fcopqr are interaction constants for the 
cubic and quartic interactions. These nonlinear interac- 
tions are generated from the Brenner potential"^^ by finite 
difference method. 



FIG. 1: Feynman diagrams for G,,(r), with the cubic and 
quartic nonlinear interactions up to the second order of h. 



By applying Wick's theorem, Gjir) can be expanded 
in terms of Hn ■ Fig. [1] gives the Feynman diagrams for 
Gjij) up to the second order of fi. For thermal expan- 
sion, the first diagram in this figure is most important, so 
we mainly consider this diagram in the following context. 
The contribution of this diagram to G'j(r) is: 



E 



j dTldT[dT'(Tln 



(6) 



where we have introduced '7],„,„(ti, r{, r") = 
fcimn<5(T{, ri)(5(T", Ti) for convenience. The contour 
ordered GF can be transformed into time domain by 
using the corresponding rules^^i r — > (i, tr), where the 
branch index cr = ± is introduced such that t = t + iea 
is on the upper (ct = +1) or lower (a = —1) contour 
branches, e is a small positive number. We find that 
GJ(i) is independent of the branch index a and time 
t. This can be understood in the sense that the system 
after thermal expansion is in thermal equilibrium state, 
and the thermal expansion effect does not depend on 
the direction of the contour. For simplicity, we use Gj 
to represent the one point GF G'^{t) in the following. 
After some derivations and using the Langreth rules (see 
appendix), we obtain the final expression for Gj: 



G, 



Imn 



inG{',n{Q)G„ 



[0], 



(7) 



where and C are the greater and retarded GF. They 
can be calculated in the eigen space without doing any 
integration as given in the appendix. For simplicity we 
have omitted the superscript 0, which indicates a har- 
monic system. To avoid confusion, we use G{t) in time 
domain and G[uj\ in frequency domain. We note that the 
information of temperature T is carried by G^. Then 
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the averaged vibrational displacement for atom j can be 
obtained from Eq. ([1]). From now on, the symbol Uj will 
denote the averaged vibrational displacement for atom 
j. This is different from its original meaning in Eq. ([5]) 
where no average is taken. The CTE at temperature T 
is calculated from: 



ih 



1 



— X — 

M Xi 



dT 



(8) 



where M is the mass of atom. Xj is the position of atom j, 
and Xj = at the common boundary of the fixed and the 
other regions. Using this equation, the CTE is obtained 
from the information of a single atom j. To obtain more 
accurate value for CTE, one needs to average aj over the 
atoms. 



III. RESULTS AND DISCUSSIONS 






Fig. [2] is the configurations for SWCNT and graphene 
in our calculation. For both systems, the smaller balls 
(blue online) on the left boundaries are fixed atoms while 
the right boundaries are free. This boundary condition 
corresponds to thermal expansion under zero pressure. 
The X axis in SWCNT is along its axial. In graphene, the 
X axis is along the in-plane horizontal direction, y in the 
vertical direction and z in the perpendicular direction. 
Periodic boundary condition is applied in the y direction 
in graphene. 
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FIG. 2: (Color online) Configurations for different systems: 
(a) Single-walled carbon nanotube; (b) Graphene sheet. The 
left boundary parts (blue online) are fixed. 

In our calculations, the system is relaxed using the 
Brenner's potential. Then the dynamical matrix and 
nonlinear force constants are obtained by finite difference 
method. The dynamical matrix is diagonalized to obtain 
the eigen frequencies and eigen vectors. Finally, the eigen 
frequencies, eigen vectors and nonlinear force constants 
are used to calculate the GF following the formula in the 
appendix. Generally, the relevant anharmonic value from 
Brenner's potential is about \kimn/kij \ Ri O.SA^^, where 
kij is one element in the dynamical matrix. 



FIG. 3: The averaged vibrational displacement u of SWCNT. 
(a), u v.s. a:: at T = 1000 K. (b). u v.s. T for the atom 
located at x = 3 A. 

We firstly study the axial thermal expansion in 
SWCNT. Fig. O (a) shows the averaged vibrational dis- 
placement in the x direction for different atoms, i.e., u{x). 
The horizontal axis is the position of each atom along x 
axis. Generally, u{x) increases with increasing x which 
reveals the fact that the thermal expansion is homoge- 
nous along X axis. On the two boundaries, there are some 
boundary effects due to the different local environments 
from atoms in the center region. Fig. [3] (b) shows the 
temperature dependence for u of one of the atoms in the 
center region. At high temperatures, w oc T is the clas- 
sical behavior. With temperature decreasing, quantum 
effect appears, and around T — 800 K, the curve deviates 
distinctly from the classical behavior. This implies that 
the Debye temperature in the SWCNT is typically about 
800 K. The quantum effect dominates the low tempera- 
ture region with T < 150 K. In this temperature region, 
u keeps a nonzero constant, resulting directly from the 
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combination effect of the thermal zero-point motion and 
nonlinear interaction in SWCNT. It means that even at 
T = K, atoms still vibrate around equilibrium posi- 
tions, and as a result of the nonlinear interaction, their 
averaged vibrational displacements is about 1% of lattice 
constant from equilibrium positions. This value is not 
small because the interaction in covalent bonding carbon 
system is very strong, so the frequencies of the optical 
phonon modes are high, thus the zero-point motion en- 
ergy is also high. We expect this temperature dependent 
curve for u can be confirmed by x-ray or other optical 
experiments. 
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FIG. 4: (Color online) CTE of SWCNT calculated from the 
thermal expansion of atoms at different location. 



I 

CD 

I 

o 



I 

CD 

I 

O 



(a) ' 

- 


1 1 




L= 15 A 




L = 25 A 




L = 30 A 

1 1 




500 1000 1500 2000 
T(K) 



FIG. 5: (Color online) Axial CTE v.s. temperature in 
SWCNT with: (a), diameter d = 4 A, and different length L; 
(b). length L = 30 A, and different diameter d. 



The axial CTE («„) of SWCNT can be calculated by 
using Uj of each atom from Eq. ([5]). Fig. S] shows the 
value of axial CTE calculated from different atoms at 
T — 1000 K. The horizontal axis is the x position for 
each atom. SWCNT in this figure is zigzag (5, 0), with 
length i = 30 A. Due to boundary effects on the left and 
right ends, values for CTE in these two regions deviate 
obviously from that in the center region. So we drop 
data on two boundaries and average over other data in 
the center region to obtain the value of CTE. 

Fig. [5] shows the temperature dependence for axial 
CTE in zigzag SWCNT of different sizes. They are sim- 
ilar to the results of armchair SWCNT. It shows that 
the length and diameter have small effects on axial CTE, 
and the value of axial CTE is positive in the whole tem- 
perature region. This is consistent with previous results 
in Ref. [l^, but different from that in Ref. UJ, where the 
axial CTE depends on diameter sensitively and is neg- 
ative at low temperature. In general, the value of CTE 
increases quickly in low temperature region, and reaches 
the classical limit at high temperatures. We point out 
that at T = 800 K, this value still changes a lot with 
changing temperature. This temperature (800 K) is a bit 
smaller than the value of Debye temperature in SWCNT 
or graphene, which is above 1000 Kj^i^i^ This result 
can be understood in the sense that the optical phonon 
modes play an important role in the calculation of the 



Debye temperaturcj^S. While in the study of the thermal 
expansion effect, the optical phonon modes are not very 
important due to their high frequency. So the critical 
temperature here is lower than Debye temperature. 

Fig. [6] is the temperature dependence for the radial 
CTE {ar) of different sized SWCNT. The radial CTE 
is obtained from the changing of diameter with chang- 
ing temperature. Compared with the axial CTE, there 
are abundant behaviors here. (1). Firstly, we can see 
Q!r < in low-temperature region. The origin of this 
phenomenon is the bending vibration in this rod-like 
SWCNT systemi^ Fig. [7] shows the position for midline 
in SWCNT along the x axis at T = 1000 K. This figure 
indicates the bending vibration directly. This bending 
vibration has a very low energy thus is very important 
at low temperatures. It leads to contraction in radial 
direction in the SWCNT. As a result, the radial CTE 
is negative in low temperature region where other high 
energy vibration modes are not excited. (2). A general 
behavior in both Fig. [6] (a) and (b) is that the thinner 
SWCNT (with higher aspect ratio length/diameter, L/d) 
has smaller value of and larger negative CTE region. 
This is mainly because the bending vibration in thin- 
ner SWCNT is more intense (see Fig. ^. For SWCNT 
with different diameters shown in Fig.[6](b), there is one 
more reason. In SWCNT, another important vibrational 
mode in the radial direction is the radial breathing mode. 
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FIG. 6; (Color online) Radial CTE of SWCNT with (a). 
d — 4: K and different length; (b). L = 30 A and different 
diameter. 



FIG. 7: (Color online) The averaged position of axial line at 
T = 1000 K in zigzag SWCNT with (a), d = 4 A and different 
length; (b). L = 30 A and different diameter. 



This mode can make the SWCNT expand in the radial 
direction. The frequency of the radial breathing mode 
is inverse proportional to diameter. So in SWCNT, the 
larger the diameter, the smaller the frequency of radial 
breathing mode. Therefore, more positive contribution 
at low temperatures will lead to larger value of a^. The 
competition between the bending vibration and the ra- 
dial breathing vibration gives a valley in all curves of ar 
in Fig. [6l (3) . The value of ar is typically one order of 
magnitude smaller than Ua- Because along the axial of 
SWCNT there are some high frequency optical phonon 
modes, which lead to strong expansion in the axial direc- 
tion. Another reason is that along the circular direction, 
SWCNT forms a closed configuration, which makes it 
more difficult to expand in the radial direction. While in 
the axial direction, the right side is open and expansion 
along this direction is easier. In Fig. [6] (a), the CTE does 
not show convergence behavior with increasing L. This 
is mainly due to very serious bending vibration in long 
SWCNT. As a result of this strong bending movement, 
it will be difficult for our method to calculate the radial 
CTE accurately. 

Now we turn to the thermal expansion in graphene. 
Graphene is easily bent due to its two-dimensional elas- 
tic sheet configuration. So in the experiment it is very 
difficult to hang a graphene sheet by just holding it on 
one side (e.g. left side as shown in Fig. [21(b)). Further- 
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FIG. 8: (Color online) CTE v.s. temperature for graphene 
sheet with various strength of interaction between substrate 
and graphene. The interaction 7 is from to 10 eV/(A^it). 



more, to study the thermal expansion effect, one can not 
suspend the graphene sheet by fixing both left and right 
sides. The only reasonable choice is to put the graphene 
sample on a substrate. In this situation, the substrate 
can support the graphene stably. But in the mean time, 
interactions between the substrate and graphene will af- 
fect the bending vibration seriously, thus influences the 
measured value for CTE in graphene sheet. In the fol- 
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lowing we focus on this substrate effect on the CTE in 
graphene sheet. 

The interaction between substrate and carbon atoms 
in graphene can be described by the on-site potential 
in the z direction^: V = (7/2)1*^ with 7 as the inter- 
action force constant in unit of eV/(A^u). is the 
vibrational displacement in the z direction. Fig. [5] is 
the CTE of graphene with 204 atoms for different val- 
ues of 7. CTE is very sensitive to 7. When there is 
no substrate interaction the bending vibration is very 
strong (which has negative effect on CTE), so the CTE 
is very small in the high temperature limit, and a large 
negative CTE region occurs at low temperatures. At 
T = 300 K, a = -6 X 10"'^ Yr^ . This value agrees 
with the very recent experimental result of —7 x 10~^ 
K~^ji However, when the substrate interaction is intro- 
duced, even a very small value of 7 = 0.001, can sig- 
nificantly enhance the value of CTE. Furthermore, the 
negative CTE region considerably shrinks. After further 
increasing interaction, the value of CTE is further in- 
creased; the negative CTE region is further reduced and 
disappears after 7 = 0.1. A saturate value for CTE is 
obtained after 7 = 1.0, where the bending vibration has 
been suppressed completely by the substrate interaction. 
We should pointout that the value of 7 = 0.001 is actu- 
ally only about 0.06% of the interaction between carbon 
atoms in the two-dimensional graphene sheet^^ which is 
more than 1.67 (N jiK^w). The practical experimental 
substrate interaction might be larger than this value. So 
the measured value of CTE will be higher than the curve 
with 7 = 0.001 in Fig. [5] If the experimental substrate in- 
teraction is larger than 7 = 0.1, no negative CTE should 
be observed. By the way, if we compare the saturated 
value at high temperatures with that in the SWCNT 
as shown in Fig. [5l we find that the value of CTE in 
graphene is much larger than SWCNT. This fact implies 
that the bending vibration in the SWCNT is also very 
important for the axial CTE. It considerably reduce the 
value of axial CTE without causing negative axial CTE 
in SWCNT. 



IV. CONCLUSION 

In conclusion, the GF method has been applied to in- 
vestigate the CTE in SWCNT and graphene. The effect 
of quantum zero-point motion and nonlinear interaction 
at T = K has been observed clearly. The axial CTE 
in SWCNT is positive in whole temperature range, while 
the radial CTE can have negative value at low tempera- 
tures. It shows that thinner SWCNT has smaller radial 
CTE, resulting from the competition between bending vi- 
bration and radial breathing vibration. Our calculation 
displays that the CTE in graphene sheet is very sensi- 
tive to the interaction between substrate and graphene, 
which will greatly enhance the CTE in graphene sheet. 
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APPENDIX A: LANGRETH RULES 

To arrive at Eq. ([7]), relations between different types of 
Green's functions are in needi^ The following Langreth 
rules^i^i are used to obtain these relations. 

For the time invariant contour integral 



(Al) 



its Fourier transform has following relations: 



A>[uj] = ^■'■[cj]C>[cj] +S>HC"[, 



(A2) 



APPENDIX B: GREEN'S FUNCTION IN EIGEN 
SPACE 

The retarded GF G^[uj] can be obtained straightfor- 
wardly, since it has the expression as^: 



GSM - [iij + i7jyi-K]-\ 



(Bl) 



where K is the N x N force constant matrix, and rj is 
a small positive number. TV is the degrees of freedom in 
the system. / is a unit matrix with the same dimension 
as K. As a result, 



G5[0] 



(B2) 



The calculation of G^(0) is more complicated. How- 
ever we find that it will be much more advantageous to 
work in the eigen space. 

(a). Firstly, the force constant matrix K is diagonal- 
ized: 



S^'KS 



(B3) 



where the unitary matrix S stores the information of all 
eigen vectors. And eigen values are in the diagonalized 
matrix f2^: 



/ 



(B4) 
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(b). The GF in real space (including both time do- 
main and frequency domain) can be obtained from the 
corresponding one in the eigen space. We take Gq[uj] as 
an example to illustrate this correspondence: 



G5M 



[{Lu + if]fi - sn^s^]^' 

S[{Lu + irifl-n^]-^S'' 



(B5) 



where we have used the relation for two matrices 
{AB)^^ = B^^A~^ and we have introduced notation 
G^iuj} = [(w + ir])^I - n^]-^ to denote the GF in the 
eigen space. Gq{uj} is a diagonalized matrix with the 
element as 



Similarly we have, 

G>{t) = S{G>{t}) Si 



(B6) 
(B7) 



This kind of correspondence is valid for other GF, and 
the relations between different GF in the real space still 
hold in the eigen space. 

(d). From Zeng's thesis^^ or P. Brouwer's note2£, the 
GF for single oscillator system with frequency ujq has 
explicit expression in time domain and frequency domain: 
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where / — /{coq) — ^pn^g_^ is the Bose distribution hmc- 
tion with 13 = l/{kBT). 

(e). Using Eq. (|B7| and ((BM)) . the one point GF Gj 
can be calculated without doing any integration: 



Gj — hmnGf^^ (O)G^, [0] 
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In the study of the coefficient of thermal expansion, 
the derivative of Gj (t) with respective to temperature T 
is needed: 
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